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Abstract 

Runaway stars are stars observed to have large peculiar velocities. Two mechanisms are thought to contribute 
to the ejection of runaway stars, both involve binarity (or higher multiplicity). In the binary supernova scenario 
a runaway star receives its velocity when its binary massive companion explodes as a supernova (SN). In the 
alternative dynamical ejection scenario, runaway stars are formed through gravitational interactions between 
stars and binaries in dense, compact clusters or cluster cores. Here we study the ejection scenario. We make use 
of extensive N-body simulations of massive clusters, as well as analytic arguments, in order to to characterize 
the expected ejection velocity distribution of runaways stars. We find the ejection velocity distribution of the 
fastest runaways {v > 80 km s^^) depends on the binary distribution in the cluster, consistent with our analytic 
toy model, whereas the distribution of lower velocity runaways appears independent of the binaries properties. 
For a realistic log constant distribution of binary separations, we find the velocity distribution to follow a simple 
power law; T{v) cx for the high velocity runaways and u"'^/^ for the low velocity ones. We calculate 

the total expected ejection rates of runaway stars from our simulated massive clusters and explore their mass 
function and their binarity. The mass function of runaway stars is biased towards high masses, and depends 
strongly on their velocity. The binarity of runaways is a decreasing function of their ejection velocity, with 
no binaries expected to be ejected with v > 150kms~^. We also find that hyper-runaways with velocities of 
hundreds of km s^^ can be dynamically ejected from stellar clusters, but only at very low rates, which cannot 
account for a significant fraction of the observed population of hyper-velocity stars in the Galactic halo. 



1. INTRODUCTION 

Runaway stars are those stars observed to have large pe- 
cuHar velocities (40 < Vr,p.r < 200kms~^ (|Blaauw 1961t 
Gies 1987: Hoogerwerf et a 1 1^0 L) and even higher ( Martini 
2006h : the specific definitions vary). A considerable frac- 



tion of the early O and B stars population are known to be 
runaway stars (~ 30 — 40 % of the O stars and 5 — 10% 
of the B stars; IS tone! [1991] and refs. therein). The velocity 
dispersion of the population of runaway s tars is much larger 
than that of the 'normal' early-type stars (lStonelll99l1) . Be- 
sides their relatively high velocities, runaway stars are also 
distinguished from the normal early-type stars by their much 
lower (< 10%) multiplicity compared with the binary frac- 
tion of normal early-type OB stars (> 50% and up to 100%; 
Garmanvetal.'"1980'; Mason e tal.ll 19981; iKobulnickv & Frveij 
2006; Kouwenhoven et al. 2007^ 



Two mechanisms are thought to contribute to the accel- 
eration of regular runaway stars, both involve binarity (or 
higher multip licity). In the binary supernova scenario (BSS; 
lBlaauwiri961l) a runaway star receives its velocity when the 
primary component of a massive binary system explodes 
as a supernova (SN). When the SN shell passes the sec- 
ondary the gravitational attraction of the primary reduces con- 
siderably, and the secondary starts to move through space 
with a velocity comparable to its original orbital veloc- 
ity. I n the dynamical ejection scenario (DES; iPoveda et al.l 
11967!) runaway stars are formed through gravitational in- 
teractions between stars in dense, compact clusters. Sim- 
ulations show that such encounters may produce runaways 
with veloc ities up t o 200km s~^, and even higher in more 
rare cases (lMikkolairi9 83; Le onard & Duncanlll990l;lLeonardl 
119911; iGualandris et alj|2004l) . These scenarios suggest that 



many of the early OB stars formed in young clusters could 
be ejected from their birth place and leave the cluster at 
high velocity. Interestingly, observations show that even very 
massive O-stars could be accelerated to become runaways 
jComeron & Pasqualill2007l) . 

In recent years, stars with extremely high peculiar ve- 
locities (hyper-velocity stars; HVSs) of a fewxlO^kms"^ 
have b ee n observed in the Galactic h alo jBrown et al.ll2005T 
l2006allR [20071: lEdelmann et al.ll2005h . HVSs are thought to 
be ejected from the Galactic center following a dynamical in- 
teraction with the massive b l ack hole (M BH ) known to ex- 
ist in the ce nter (Hills 1988; Hansen & Mlosavlievic '20031: 
lYu & Tremaine 2003; Perets et al,, 2007; Perets. 2009a) . Nev- 
ertheless, it was suggested that the BSS and/or the DES could 
eject, under some conditions, stars with extreme velocities 
(termed hyper-runaways), comparable with those of observed 
HVSs, possibly explaining the origin of some of these HVSs 
([Gyaramadze et al. 2009). 

Several studi es have explored the DES (see 

iHoogerwerfetaP [2000l) using N-body simulations. Most 
of these have focused on single encounters between a 
single/binary star and another binary star, and found the 
velocities of the eject ed stars in such encounter s (e.g . 
iGvaramadze et al.l l2009h . [Leonard & Duncan! (119901 119881) 
explored the dynamical properties of runaway stars and 
their distribution in a stellar cluster environment, and not 
in single encounters. Their studies, however, were limited 
to a relatively small number of simulation of very small 
clusters (^ 30 stars; in addition to some hybrid N-body - 
Monte-Carlo simulations of larger clusters of a few hundreds 
of stars), and he nce to very small sample of ejec ted runaway 
stars. Recently, iFujii & Portegies ZwartI (1201 1!) studied the 
ejection of stars from large clusters, however, they considered 
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only dynamically formed binaries, and did not include 
primordial binaries in their simulations. They also run only a 
limited number of simulations, which do not provide enough 
statistics of the high velocity tail distribution of runaways. 

Here we develop an analytic understanding of velocity dis- 
tribution of the fastest runways in the DES. We compare 
it with a large sample of simulated runaway stars produced 
in extensive N-body simulations of hundreds of large stellar 
clusters. These provide us, for the first time, large enough 
database of runaway stars with high velocities, which could 
be analyzed statistically. Using our N-body simulations we 
characterize the ejection rates, velocity distribution, and bi- 
narity of dynamically ejected runaway stars. We find cases 
of stars dynamically ejected at extreme ejection velocities of 
a few hundreds of kms^^, however these are relatively rare 
cases and are not likely to explain the vast majority of hyper- 
velocity stars observed and inferred to exist in the Galactic 
halo. 

2. ANALYTICAL ESTIMATES 

We begin by exploring analytically the dynamical ejection 
scenario for runaway and hyper-runaway stars. We consider 
an interaction of a binary of mass Mb = Mi + M2 and a 
single star Mi,. Large accelerations of the single star are as- 
sumed to occur when it passes one of the binary components 
within the semi-major axis, a, of the binary. Assuming that 
the scattering is dominated by the gravitational focusing, the 
cross section of the interaction is 



a{a) 



2TTGMBa 



(1) 



where Vc is the characteristic stellar velocity in the cluster 

The energy exchange between a hard binary and star is 
comparable to the binary orbital energy. Hence, the energy 
transfer to the star is of the order 



GMbM, 



(2) 



and the star acquires a large velocity kick, Wkick of the order 
of the orbital velocity of the binary. We are interested in the 
cases where AEi, is much larger that the star's energy before 
the interaction. The velocity of the ejected star Wkick — V AE 
is therefore estimated as 



^ 2GMb 



(3) 



with a = 1/2. Simulations of binary-single star encounters 
su ggests a somewhat steep er slope of a ~ 3/5 (see figure 4 
in 



Gvaramadze 



ewhat steep e 
etani2009h . 



The differential cross section per unit volume of the binary- 
single star interaction with semi-major axis in the interval 

(a, a + So) is 



dRja 
da 



n{a)NBni,a{a)vc 



(4) 



where ub is the number of binaries per unit volume, n^, is 
number density of stars and n{a) is the differential distribu- 
tion of semi-major axis of the binaries, which is normalized 
to unity: 



dan{a) = 1 



(5) 



Let us consider a log constant distribution of the binaries 
semi-major axis (so called Opik's law, with n{a) oc a~^\ 
(3 = 1), representing an empirical distribution of massive bi- 
naries (e.g. Kobulnickv & Frver 2006) and, for testing pur- 
poses, an uniform distribution n{a) = const; /? = 0. Due 
to the gravitational focusing, the cross section of encounters 
with massive binaries is much higher than with low mass 
stars. In addition, the ejection velocities from encounters 
with massive stars are higher Observationally, massive stars 
also have a much higher binarity fraction, as well as separa- 
tion distribution which is biases toward closer binaries (lower 
mas s stars seem to have a log normal distribution of peri- 
ods; iDuquennoy & Mayoilll99ll) . For all these reasons, the 
majority of runaway stars are ejected through interactions of 
massive stars with massive binaries. For simplicity, we ne- 
glect dependence of dR{a) on the specific mass function of 
binaries in Eq. and assume some fiducial effective typical 
mass for the stars involved in the dynamical encounter. In- 
terestingly, our numerical simulation results suggest that this 
may be justified even for a range of stellar masses, however, 
we leave the more detailed analytic study of this dependence 
for further investigation. 

Eq. (|4]i can be directly transformed to the production rate of 
high velocity stars per unit volume and velocity interval dvcsc- 

(6) 

Making use of Eqs. ([T]i and ^ we can find the expected ve- 
locity distribution of the ejected stars, for some appropriate 
choice of the binaries separation distribution. For an Opik's 
distribution of the semi-major axis 7i{a) oc we expect the 

velocities of ejected stars to be distributed as i?(uesc) Ucsc^"^ 
{a ~ 3/5; (3 = 1), while for n{a) = const, we get 

R{vcsc) oc fcsc'^^'^ (a = 3/5; /3 = 0). This approach can be 
naturally extended to any desired distribution of binary sepa- 
rations. 

3. N-BODY SIMULATION MODELS 

We have studied several different numerical models of star 
clusters in order to characterize the properties of runaway 
stars. We also make use of simplified cluster models to ver- 
ify our analytic calculations of the the velocity distribution 
of ejected runways. Table [T] summarizes basic physical pa- 
rameters of the models which were integrated by means of 
NBODY6 code ( Aarseth 1999). In all cases, distribution of 
the binary semi-major axis was truncated outside the interval 
(0.05 AU, 50 AU); initial eccentricities were set to zero. In 
the first models we studied a highly simplified stellar systems 
in which only two masses were used rather than a continu- 
ous mass function. All binaries in these models are made of 
identical 10 Mq components and the single stars are identical 
1 Mq stars; in addition, all stars are treated as point-masses in 
these models. Later we also compared these to simulations re- 
sults of similar clusters with no primordial binaries, to verify 
that single-single encounters can only lead to ejections with 
Vcsc ^ 20kms~^. During the cluster evolution, dynamically 
formed binaries could eject stars at higher velocities, but the 
overall fraction of runaways was only a small fraction than in 
the case of clusters including primordial binary population. 

Though in this study we focus on the simple cluster sim- 
ulations, we also run simulations of somewhat more realistic 
clusters, in which the full ranges of Salpeter initial mass func- 
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tion (IMF) for the stars is considered. In the latter models, all 
massive (M > 4 Af©) stars reside in primordial binaries, and 
we set a lower limit for the primary and the secondary star 
to 4 Mq and 1 Mq, respectively. We use pairing algorithm 
which prefers similar masses of the two components which is 
in accord with observations. Furthermore, binary stars in real 
clusters may physically collide. Hence, we enable the possi- 
bility of stellar collisions for these models in order to increase 
their realism; stars are merged if they pass to each other at a 
distance smaller than the sum of their radii. 

The process under the consideration, i.e. acceleration of 
stars to velocities > 60kms^^ is quite rare. Therefore, 
we have integrated several hundreds of different realizations 
(-^run) of each model in order to obtain statistically relevant 
results. 

4. RESULTS 

4.1. Velocity distribution 

Figure [T] shows the velocity distribution of escaping stars 
(found at least 5 pc away from the host cluster at the end 
of the simulations). The figure shows the velocity distribu- 
tion compared with the predicted distribution at the high ve- 
locity regime (when taking the overall normalization to be a 
free parameter). Although the analytic derivation does not 
account for the mass function, we find that the velocity dis- 
tribution of ejected stars in models 009 and 014 are also con- 
sistent with simple analytic formulation, and the slope is de- 
termined mainly by the distribution of the semi-major axis 
of the binaries (see fig. [T]i- The velocity distribution at the 
lower velocity regime v^sc ~ 20 — 150 km s^^ is qualita- 
tively different, and appears to be independent of the distri- 
bution of the binary separation. Comparison with a model 
where no binaries exist show that all the bona-fide runaway 
stars arise from the existence of binaries in the clusters. Al- 
though not the focus of this paper it is interesting to note in 
passing that stars with velocities in the intermediate velocity 
range 20 — 150 km s^^ can not be explained by single-single 
encounters alone; however, they also do not follow the sim- 
ple analytic approach described before, but rather appear to 
follow a shallower distribution which is independent of the 
binary separation distribution. This may be cons istent with 
the recent results of lFuiii & Portegies Zwarli (1201 ll) who stud- 
ied runways from massive clusters in this range of veloci- 
ties. These may arise from the dynamical interactions with a 
small number or even single very massive binaries in the clus- 
ter (bullies), which dominate the ejection rate and are later 
on ejected themselves from the clusters. The properties of 
runaways from such interactions are independent of the over- 
all binary population characteristics. Our findings in Fig. |4] 
showing ejection of very massive runway binaries is consis- 
tent with this picture. Note, however, that models 016 and 
017 did not include massive binaries, which can play the role 
of 'bullies'. Unfortunately, we did not keep data on specific 
interactions leading to runaway ejections; we therefore can 
not directly confirm or refute this scenario for the low veloc- 
ity regime. We defer such simulations and analysis to future 
work. 

4.2. Binarity of runaway stars 

The binarity of runway stars is an important signature of 
the dynamical ejection scenario. Previous studies of dynam- 
ically ejected runaway stars suggested their binary fraction 
to be low relative to the binary fraction of their parent stars 
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Fig. 1 . — The velocity distribution of runaway stars ejected from the simu- 
lated clusters with velocities larger than 20kms^^. The runaways velocity 
distribution for cluster of different binary separation distribution are shown; 
Dash-dotted and solid stair lines are for a log-constant SMA distribution 
(models 009 and 016, respectively) and dash and dotted lines are for uniform 
SMA distribution (models 014 and 017, respectively). See Table 1 for initial 
conditions of all simulated clusters. The right segment of the straight lines 
(solid brown and red) show the toy-model predictions for the runaway veloc- 
ity distributions, normalized to fit simulated data. Left segment of the straight 
lines (solid blue) show power law distributions with n (x v^^'^. Note that 
although the predicted slopes provide the correct trends for the models, the 
best fitting (not shown) power law slopes for the same high velocity regimes, 
V > 150kms~^, are —2.5 and —2.9 for models 009 and 016, respectively 
(compared with the predicted 7 = —2.7); and for the other models, with 
velocity regime v > SOkms"'^, we get —3.4 and —3.3, for models 014 and 
017. respectively (compared with the predicted 7 = —4.3) 

in the cluster ([Leonard & Duncanlll990l) . We find a similar 
trend; Fig. |2] shows the velocity distribution of single and 
binary runaways in models 009 and 014 of our simulations 
(for which a realistic binary period distribution was used), 
and the binary fraction of runaways as a function of their ve- 
locity. The simulations results are consistent with theoretical 
arguments; in a binary single encounter the binary receives 
a smaller fracti on of the k inetic energy, il/,t/(A/B + M*) 
(fPortegies Zwar t et alj|2010l) , it is therefore typically ejected 
at considerably lower velocities. Note that encounters where 
the single star is much more massive than the binary com- 
ponents likely lead to an exchange of the massive compo- 
nents with one of the lower mass binary components, there- 
fore even such encounters would eventually lead to the bi- 
nary being the most massive component in the encounter, fur- 
ther contributing to its low ejection velocity, compared to the 
ejection of single stars. Binary-binary encounters may eject a 
binary to higher velocities, however, these complex encoun- 
ters may easily disrupt one of the binaries and/or form higher 
multiplicity systems. A hard binary can be ejected similar 
to the ejection of single star by a softer binary, however its 
typical ejection velocity would be comparable to the orbital 
velocity of the wide binary, again producing a bias toward 
low velocity ejections. The overall binary fraction of run- 
aways decreases significantly with higher ejection velocities; 
and effectively no binary was ejected at velocities higher than 
- 150 (300) kms-i in the 009 (016) models (compare with 
the ejection velocities of single stars). Hyper-runaway bina- 
ries are t herefore not likely to be produced through the DES 
(see also lPeretil2009bl; iBrown et al.ll2010l for discussions on 
this issue). 



4 Perets & Subr 



model 


A^run 


rh (pc) 




S 


Mc (Mq) 


IMF 


n(a) 


A^bin 


Mp(M0) 


M, (Mq) 


009 


200 


0.2 


2.7 


0.25 


5000 


Salpeter 


a-i 


73 


> 4 


> 1 


014 


500 


0.2 


4.4 


0.25 


5000 


Salpeter 


const. 


« 73 


> 4 


> 1 


016 


200 


0.4 


35.5 





10000 


200 X 10 Mq + 8000 Mq 


a-i 


100 


10 


10 


017 


200 


0.4 


35.5 





10000 


200 X 10 Mq +8000Mq 


const. 


100 


10 


10 



TABLE 1 

Parameters of the models: initial half-mass radius r^, index of the initial mass segregation S (see lSubr et al.l|200^ , total mass of the cluster M^, distribution 
of the binary semi-major axis n(a), initial number of binaries /j,;^ and mass of the primary A/p and secondary A/g star. Salpeter initial mass function has 

power-law profile oc M^^ '^^ within the interval (0.2 Mq, 80 Mq). In the models 009 and 014, physical collisions of stars are allowed (i.e. the mass function 
is not constant). Initially, the binaries have semi-major axes in the range (0.05 AU, 50 AU) and zero eccentricities. Profiles of initially mass segregated models 
cannot be expressed analytically but, in general, their density increases towards the center 
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Fig. 2. — The binary properties of runaway stars. The binary fraction of 
runway stars, as a function of the ejection velocity in models 009 and 016. 
Inset: The cumulative velocity distribution of runaway binary stars. 

4.3. The mass function and mass-velocity relation for 
runaway stars 

In Fig. |3]we show the runaway fraction of stars. As can be 
seen, OB runaways are relatively more abundant than lower 
mass stars, with the O-star runaways fraction 2-3 times larger 
than that o f the B-stars, consistent with observations (see 
IStonelll991l andrefs. therein). 

The total fraction of runaways we find are comparable to, 
but sys tematically lower (~ 1/2) than those reported by Stone 
(1199 Ih : this may result from the limited time of the simula- 
tion; the lifetime of typical O-stars could be 2-3 times longer 
than the simulation time (which is comparable to the life- 
time of the most massive stars in our simulations). Obviously, 
many other simplified assumptions we use affect our theoret- 
ical results, and may contribute to the difference. 

In Fig. |4]we present the mass-velocity distribution of the 
ejected stars. We see again the trend of more massive stars 
having higher runaway fraction (as reflected by the large frac- 
tion of massive runaways compared with their fraction in the 
IMF); we note that even the most massive stars are ejected 
as runaways. We also find that the velocity distribution of 
more massive runaways tends towards higher velocities than 
the lower mass stars, i.e. the fraction of more massive stars 
increases with ejection velocity. In fact, we find that ^ 20 
(10) % of the runaway O-stars with v > 20 (100) km s^^ 
are more massive than 40 A/©, and ~ 5 (3) % of the O-stars 
are with u > 20 (100) kms^ are more massive than 80 M© 
(consistent with finding of very massive runaway stars; e.g. 
iHoogerwerf et al.ll2001h . The latter are runaway merger prod- 
ucts of two or more less massive stars, since the adopted IMF 
of the clusters upper mass cut-off was 80 M©. 
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Fig. 3. — The fraction of runaways vs. mass with given (lower limit) ve- 
locities for model 009. All stars below 4 Mq are initially single stars, and 
their runway fraction might be non-realistic. For OB stars m > 4 Mq the 
runway fraction increases with mass. This trend is consistent with the ob- 
served runaway fractions of O stars to be 3 — 6 times larger than the runaway 
B-stars fraction (see Stone ( 1991) and refs. therein). Note that even the most 
massive stars are ejected as runaways; some of which are merger products, 
with masses extending beyond the initial mass distribution (not shown). 
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Fig. 4. — The mass function of runaway stars in mod 009 as a function 
of lower limit velocity. The distribution tends to be bi-modal, with runways 
in the range 1 — 4 Mq typically being least represented in the population 
(compare with the shown IMF of all stars in the initial simulated clusters). 

4.4. Hyper-runaways and hyper-velocity stars 

As can be seen for both the two-mass models and the con- 
tinuous mass ones, the high velocity tail of the distribution 
follows a steep power law. The ejection or the fastest run- 
aways requires the close encounters with binaries, the rate 
of which are dominated by the most dense clusters (which 
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are also the most massive ones, typically). Therefore, in or- 
der to provide a basic estimate of the ejection of the fastest 
runaways we consider here only the most massive clusters in 
the Galaxy. The average number of hyper-runaways ejected 
with Uosc > 300 (> 450) kms^^, in our most realistic clus- 
ter, model 009, is ^ 0.2 (0.02) stars per cluster Currently, 
Nc ~ 10 — 20 young (< = 10 Myrs) massive clus- 
ters {M > IO^Mq) with cores comparable to our simu - 
lated clusters exist in the Galaxy jMurrav & Rahmanll2010l) . 
We therefore expect the ejection rate of hyper-runways to be 
- 0.2 X iV, X = 0.2 - 0.4 Myr-i. For B-stars (of 
masses ~ 3 — 1(3 M0 such as observed among the hyper- 
velocity stars in the Galactic halo), the rate is ^ 2.5 times 
lower Over 100 Myrs (comparable to the propagation time 
for the observed HVSs in the Galactic halo), we therefore ex- 
pect to have of order 10-20 hyper-runaways in the galaxy, but 
only one or two such stars with Uosc > 400. Currently ^ 20 
HVSs (wcsc > 400; though this is generally a good definition, 
note that a more accu rate definition also depends on distance 
from the Galaxy, see iPerets et all 120091) have been observed 
in the Galactic halo, from which a total of ^ 10 B-stars (of 
3 — 4 M(7)) are inferred to exist in the Galaxy jBrown et al.l 
I2006ai 12009 ; Brown 2()1JJ, and references therein), and a 
few hundred lower v elocity HVSs (so c alled bound HVSs; 
275 < i^esc < 400 iBrown et al.l l2007h . may exist in the 
Galaxy (iPerets et al.ll2009h . Comparing the predicted and ob- 
served numbers of HVSs, hyper-runaways from this scenario 
may contribute at most 1 — 2 % of the HVSs population and 
are unlikely to have produced any of the currently observed 
HVSs in the Galactic halo (or at most one); though they may 
contribute a small fraction of the observed bound HVSs. 



distribution of the fastest runaway stars. We then used ex- 
tensive N-body simulations of simplified clusters composed 
of only two type of stars with masses 1 and 10 Mq to char- 
acterize the velocity distribution of runway binaries; later we 
also simulated clusters with realistic continuous Salpeter IMF. 
We found the trends in velocity distribution to depend on the 
properties of binaries in the cluster and be generally consis- 
tent with the analytical toy models of multi-interactions with 
multiple binaries at the high velocity regime v > 80kms~^ 
(or V > 150 km s^^; depending on the model). At the lower 
velocity regime, ejection velocities are less dependent or pos- 
sibly even independent of the binary properties, suggesting a 
diffe rent channel for runaway eject ions dominate this regime 
(e.g. iFuiii & Portegies Zwarfcol lb . 

We characterized the velocity distribution of runways from 
the continuous IMF clusters, and discussed their mass func- 
tion and binarity. We find that the runaway fraction of stars in- 
creases with mass, consistent with observations; and that very 
massive stars formed through collisions in the cluster could 
be ejected as runaways. The binarity of runaways stars de- 
crease with increasing velocity, and is generally lower than 
the binary fraction of their birth cluster In particular, we 
find the maximal velocity of binary runaways is limited to 
< 200 km s^^. We also find the mass function of runaways 
to be velocity dependent. Although runaways are much more 
frequent amongst the massive stellar population > 4 Mq, we 
find that a large population of low mass runaways should also 
exist. The ejection rate of hyper-runaways, with velocities 
> 300 km s^^ appear to be too low to explain a significant 
fraction of the observed HVSs in the Galactic halo, and could 
at most explain a small fraction of the observed bound HVSs. 



5. SUMMARY 

In this paper we explored the dynamical ejection of run- 
away stars from the cores of massive stellar clusters. We pre- 
sented a simple analytical toy model to explain the velocity 
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